function bgk_modes

%fname="kx010jv560_dt1e-2.bgk";  fbase="bgk_kz035_E"

fname="jv0800_dt1.dfu.0001.00";  fbase="jv0800_Dv_E"

%jv1200_dt1_free.dfu.0003.00
%jv1200_dt2.dfu.0001.00
%jv1200_dt2_free.dfu.0003.00
%jv2400_dt2.dfu.0001.00
%jv2400_dt2_free.dfu.0003.00


%---------------------------------------------------------

% system size

jz=64; jv=800;

% physics paramters

vmax=8.00;  
kz=0.35;
#vPhase=3.4;
vPhase = 0;

iunit = i;


%---------------------------------------------------------

  nz=1;
  dkz = kz/nz;
  dz = 2*pi*nz/(jz*kz);
  z  = (0:jz-1)*dz;
  Lz = 2*pi/kz;

  dv = 2*vmax/jv;
  dp = 2*pi/(2*vmax);
  p  = fftshift([-jv/2:jv/2-1],2)*dp;  %p(ind)=0;
  kz = fftshift([-jz/2:jz/2-1],2)*dkz;

  v  = (0:jv-1)*dv - vmax - vPhase;
  f0 = ( 1/sqrt(2*pi) )*exp( - 0.5*(v + vPhase).^2 );

  fbg = meshgrid(f0, 1:jz)';  %'
  pp  = meshgrid(p, 1:jz)';   %'

  

%---------------------------------------------------------

  fid = fopen(fname, 'rb');
  dfu = fread(fid, jv*jz, 'double');
  fclose(fid);

  dfu = reshape(dfu, jv,jz);
  %f = dfu + fbg;
  f = dfu;
  %imagesc(flipud(f));   axis('equal', 'tics', 'off');
semilogy(v, abs(f(:,34)));   
axis([-8,8,1.e-16,1.e-1]); xlabel("kz"); ylabel("f_p");

 

  return

  %-- analytical continuation --

  %ff = ifftn(f);
  %ff = ff .* exp(-pp);
  %f  = fftn(f);

  %imagesc(flipud(log(abs(fftshift(ff)))));
  %semilogy(-jv/2:jv/2-1, fftshift(abs(ff(:,1)),1)); 
  %axis([-400,400,1.e-20:1]); xlabel("kz"); ylabel("f_p");
  %return

  %-- derivative --

  ff = fft(f,jv,1);
  fD = iunit*ff.*pp;
  Df = ifft(fD,jv,1);


  %-- modes --

  fDf = ifft(Df,jz,2);
  ff  = ifft(f,jz,2);


  %-- electric field --

  dens = sum(f,1)*dv;
  fdens = ifft(dens);
  fE = - iunit*fdens./kz; fE(1)=0;

  

% plot(v,Df(:,1));  axis([-6,6, -0.04, 0.04]); xlabel("v"); ylabel("df/dv (z=0)");
% plot(v,Df(:,1));  axis([1,4, -0.04, 0.04]); xlabel("v"); ylabel("df/dv (z=0)");
% plot(v,Df(:,1));  axis([2,3, -0.04, 0.04]); xlabel("v"); ylabel("df/dv (z=0)");


  %plot(z, dens); axis([0,Lz]); xlabel("z"); ylabel("rho");
  %semilogy(0:31, abs(fE(1:32)), 'o-' ); axis([0,25, 1.e-10, 1]); xlabel("m"); ylabel("|E_{zk}|");


  %return

%-----------------------------


  %imagesc(flipud(dfu));  axis('equal', 'tics', 'off');
  %imagesc(flipud(abs(fdfu)));  axis('equal', 'tics', 'off');
  %imagesc(flipud(Df)); axis('equal', 'tics', 'off');

  %plot(v, abs(ffbg(:,1)) );
 
  %semilogy(v, dfu_avg+f0'); axis([-6.25,6.25,1.e-4,1]); xlabel("v"); ylabel("(f0+df)_{avg}");%'
  %plot(v, 100*dfu_avg, v, f0); axis([-6.25,6.25]); xlabel("v"); ylabel("f0,   100 df_{avg}");%

  %plot(v,real(fdfu_sum), v,real(ffbg_sum)); 
  %xlabel("v"); ylabel("real(sum(df_k)),  real(f0_k)"); axis([-6.25,6.25]);

  %plot(v,imag(fdfu_sum), v,imag(ffbg_sum)); 
  %xlabel("v"); ylabel("imag(sum(df_k)),  imagbgk(f0_k)"); axis([-6.25,6.25]);


  %plot(v,real(fdfu_sum), v, real(fdfu(:,1:4))); 
  %xlabel("v"); ylabel("real(sum(df_k)), real(dk_{0,1,2,3})"); axis([-6.25,6.25]);

  %plot(v, imag(fdfu_sum), v, imag(fdfu(:,1:4))); 
  %xlabel("v"); ylabel("imag(sum(df_k)),  imag(dk_{0,1,2,3})"); axis([-6.25,6.25]);

 
 
 %-- write spectrum of the electric field --

  fname = [fbase, '.txt'];
  fid = fopen(fname, 'wt');
  fprintf(fid,'# 1.m   2.real  3.imag\n\n ')
  for iz=1:jz/2
     fprintf(fid, "%4d  %20.12e  %20.12e\n", ...
  	     iz-1, real(fE(iz)),  imag(fE(iz)) );
  end
  fclose(fid);
  
  return



 %-- write spectra of bgk deviation --

  total = ff;
  %total = fDf;

  for m = 0:31
     fn = num2str(m,'%02d');
     fname = [fbase, '_m', fn, '.txt'];
     fid = fopen(fname, 'wt');

     %fprintf(fid,'# 1.v  2.real  3.imag\n\n ')
     for iv=1:jv
        fprintf(fid, "%14.6e  %20.12e  %20.12e\n", ...
	       v(iv), real(total(iv,m+1)),  imag(total(iv,m+1)) );
     end
     fclose(fid);
  end    

end
